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SUMMARY 

We  present  a  new  method  for  stochastic  shape  optimisation  of  engineering  structures.  The  method  generalises 
an  existing  deterministic  scheme,  in  which  the  structure  is  represented  by  a  level-set  method,  and  evolves  by  steepest 
descent  of  the  objective  function.  In  non-convex  optimisation  problems,  the  deterministic  algorithm  can  get  trapped 
in  local  optima:  the  stochastic  generalisation  enables  sampling  of  multiple  local  optima,  which  aids  the  search  for 
the  globally- optimal  structure.  The  method  is  demonstrated  for  several  simple  geometrical  problems,  and  a  proof-of- 
principle  calculation  is  shown  for  a  simple  engineering  structure. 


I.  INTRODUCTION 


Topology  optimisation  has  been  demonstrated  to  provide  substantial  performance  improvements  and/or  weight 
savings  in  a  wide  range  of  engineering  design  problems  mm-  However,  many  of  the  relevant  applications  involve 
non-convex  optimisation  problems  with  multiple  locally-optimal  designs:  these  might  be  associated  with  different 
topologies  or  with  numerical  aspects  of  the  (discretised)  computational  problem.  Examples  of  non-convex  design 
spaces  are  stress  constrained  optimization  |5]  and  coupled  aeroelastic  wing  optimization  [5j. 

Non-convexity  and  multiple  optima  present  interesting  dilemmas  to  engineers.  When  a  local  optimum  solution  is 
significantly  worse  than  the  global  optimum,  the  optimisation  scheme  fails.  However,  locally-optimal  solutions  can 
also  provide  multiple  design  ideas  to  engineers,  offering  a  range  of  possible  solutions  that  an  engineer  may  wish  to 
consider,  based  on  practical  design  requirements  such  as  ease  of  manufacturing.  This  particularly  true  when  there  are 
several  designs  of  similar  objective  function  and  constraint  values,  an  example  of  which  can  be  found  in  0. 

A  common  approach  to  topology  optimisation  is  to  employ  gradient-based  nonlinear  programming,  which  can  be 
applied  to  typical  engineering  design  problems  with  104  —  106  design  variables  mm-  In  this  case,  iteration  of  the 
optimiser  leads  to  one  local  optimum  solution.  Alternative  approaches  such  as  evolutionary  algorithms,  particle 
swarm  optimisation  and  simulated  annealing  are  capable  of  searching  for  multiple  potential  solutions  and  they  have 
been  applied  to  topology  optimisation  JT2]  !2H  Si  •  However  the  success  of  such  methods  has  been  limited,  partly 
because  they  do  not  typically  take  advantage  of  information  about  the  gradient  of  the  objective  function.  Interested 
readers  are  referred  to  a  critical  review  [T5]  which  presents  an  example:  Topology  optimisation  was  applied  to  a 
small  problem  with  just  144  variables,  which  was  solved  using  the  non-gradient  method  of  differential  evolution.  This 
required  15,730  function  evaluations.  In  contrast,  a  gradient-based  topology  optimisation  method  Solid  Isotropic 
Material  with  Penalisation  (SIMP)  -  converged  to  a  slightly  superior  solution  after  just  60  function  evaluations.  This 
motivates  the  research  question  of  how  to  explore  non-convex  engineering  design  spaces  effectively  and  efficiently. 

One  approach  which  can  be  successful  in  non-convex  optimisation  problems  is  to  start  with  a  deterministic  optimi¬ 
sation  procedure  and  to  add  a  stochastic  component,  so  that  the  objective  function  can  both  increase  and  decrease  as 
the  algorithm  runs.  In  this  report,  we  introduce  a  stochastic  optimisation  procedure  which  aims  to  generate  designs 
according  to  a  prescribed  distribution,  analogous  to  the  Gibbs-Boltzmann  distribution  in  physics.  The  method  is 
built  on  the  deterministic  gradient-based  topology  optimisation  method  of  Dunning  and  Kim  mm  and  that  method 
is  recovered  in  the  zero-noise  limit  of  the  stochastic  process  -  this  means  that  the  stochastic  method  should  perform 
at  least  as  effectively  as  the  deterministic  one.  Moreoever,  since  the  model  is  based  on  the  Boltzmann  distribution, 
we  expect  that  it  can  be  combined  with  parallel-tempering  methods  m  mi  m  which  offer  a  systematic  approach  for 
exploring  a  range  of  near-optimal  structures. 

The  new  method  is  based  on  a  recently-introduced  deterministic  method  for  level-set  topology  optimization 
(lsto)  mm-  Unlike  traditional  topology  optimisation  (e.g.  SIMP)  in  which  the  design  variables  (for  example  n) 
indicate  whether  each  element  of  the  design  domain  should  exist  ( n  =  1)  or  be  absent  {n  =  0),  the  level  set  method 
employs  an  implicit  functional  representation  (j>  that  directly  represents  the  domain  boundaries.  A  key  advantage  of 
this  method  is  that  topological  changes  in  the  shape  of  the  object  (for  example  the  removal  of  holes)  are  associated 
with  singularities  in  the  behaviour  of  the  boundary  of  the  object,  but  do  not  involve  any  singularities  in  </.  See  m 
for  a  description  of  the  level-set  method  and  its  key  advantages.  The  effect  of  the  level-set  representation  is  that 
topology  optimisation  can  be  reformulated  as  an  extended  shape  optimisation,  where  boundary  shapes  are  optimised 
and  the  number  of  boundaries  can  change.  In  this  paper,  we  will  focus  on  the  shape  optimisation  aspect  of  the 
algorithm,  where  the  level  set  method  moves  the  boundaries  using  an  advection  equation.  Within  this  scheme,  the 
LSTO  method  corresponds  to  steepest  descent  of  the  objective  function,  which  ensures  robust  convergence  of  the 
method  to  a  locally-optimal  design.  We  will  show  how  this  method  can  be  extended  to  a  stochastic  method  that  can 
explore  non-convex  design  landscapes. 

In  this  paper,  we  present  several  new  results.  In  Sec.  [TIJ  we  review  the  method  of  [7],  and  we  introduce  some 
simplications  to  that  method,  which  clarify  the  relationship  between  that  method  and  steepest  descent  optimisation  of 


the  objective  function.  In  Sec.  Ill  we  introduce  a  stochastic  process  that  explores  a  range  of  structures,  parameterised 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


4 


Discretiscd  structure 


Discretiscd  boundary 


FIG.  1.  The  level-set  domain:  (Left)  A  two-dimensional  square  domain  fid  is  discretised  using  a  uniform  square  grid.  A 
function  (f>  defined  on  the  nodes  of  the  grid  indicates  whether  each  node  is  inside  (c f>  <  0)  or  outside  (4>  >  0)  a  (circular) 
structure  fl.  (Middle)  The  function  (f)  is  given  by  the  signed  distance  to  the  nearest  point  on  the  boundary  of  f 1.  With  this 
choice  |V<(>|  =  1.  (Right)  A  discretised  representation  of  the  boundary  of  f 1  is  obtained  by  defining  boundary  points,  which  are 
located  either  on  nodes  (if  4>i  =  0)  or  on  edges  between  nodes  (if  this  edge  connects  two  nodes  with  opposite  signs  for  <(>).  The 
resulting  set  of  boundary  points  provide  a  piecewise  linear  approximation  for  of  the  boundary  of  f 2. 


by  a  noise  strength  T  that  is  analogous  to  the  temperature  in  statistical  physics.  Sec.  |IV|  includes  several  examples  of 
the  application  of  the  stochastic  method,  including  matching  of  a  shape  to  a  fixed  design,  and  compliance  minimisation 
of  a  simple  two-dimensional  engineering  structure,  subject  to  a  constraint  on  its  area.  We  discuss  the  nature  of  the 
shapes/structures  explored  by  the  stochastic  method,  and  discuss  the  potential  use  of  the  method  for  practical 
optimiation  of  engineering  structures.  Sec.  |V|  summarises  our  conclusions. 


II.  THE  LEVEL  SET  TOPOLOGY  OPTIMISATION  METHOD 


The  optimisation  problem  considered  here  is  the  determinisation  of  a  structural  domain  Q  that  minimises  an 
objective  function  F(Q),  subject  to  some  constraints.  This  section  introduces  the  computational  method  used  here, 
which  is  closely  based  on  the  method  of  Dunning  and  Kim  [7] ,  see  also  m-  We  refer  to  that  method  as  the  Level-Set 
Topology  Optimisation  (LSTO)  method.  The  implementation  described  here  is  a  slightly  simplified  version  of  the 
original  LSTO  method:  a  C++  code  for  the  method  discussed  here  will  be  made  available  on  publication  of  these 
results.  We  give  a  brief  and  informal  description  of  this  (deterministic)  optimisation  algorithm,  to  set  the  scene  for 
the  stochastic  method  that  we  will  describe  in  Sec.  m  Rigorous  discussions  of  the  properties  of  the  level-set  method 
and  of  the  mathematical  results  underlying  this  algorithm  can  be  found  elsewhere  |15j. 

We  introduce  an  objective  function  F ,  a  constraint  function  G,  and  a  design  domain  f C  Rd,  where  d  is  the  spatial 
dimensionality  (In  this  work  we  take  d  —  2  although  generalisation  to  higher  dimension  is  possible.)  The  aim  of  our 
optimisation  is  find  a  structure  D  C  fid  such  that  F(fl)  is  minimised,  subject  to  the  constraint 

G(fi)  >  G*.  (1) 

Extension  to  multiple  constraints  or  equality  constraints  of  the  form  G(fl)  =  G*  is  straightforward  [7,  [TO]  but  we 
consider  a  single  inequality  constraint  here,  for  simplicity.  The  LSTO  method  prescribes  a  time-evolution  for  0  that 
converges  to  a  (local)  optimum  of  F,  which  satisfies  the  constraint. 


A.  Evolution  in  continuous  space  and  time 


The  structure  SI  is  defined  in  terms  of  a  real- valued  function  cf> :  fid  — >  R.  That  is,  define  0  =  {x  £  Dd  :  <f>( x)  >  0} 
as  the  part  of  the  domain  for  which  (j>  >  0.  Throughout  this  work,  bold  vectors  such  as  x  indicate  vectors  in  Rd.  The 
boundary  of  fl  is  denoted  by  T  and  is  defined  as  the  zero  level-set  of  4>,  that  is  T  =  {x  £  fid  :  <f>{x)  =  0}. 

The  boundary  T  is  made  up  from  one  or  more  closed  curves,  so  we  index  the  points  in  T  by  an  internal  co-ordinate 
u  >  0,  such  that  X(u)  £  T  is  a  point  on  the  boundary  of  f 1.  Also  let  n(u)  be  the  (outward)  normal  vector  to  the 
boundary  at  the  point  X(u),  and  define  a  function  £'  such  that  /“2  d£(u)  =  /“2  £'(u)du  is  the  length  of  the  boundary 
between  the  two  points  X(ui)  and  X(u2)- 
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Several  level-set  topology  optimisation  approaches  Hi  use  a  steepest  descent  strategy  for  the  minimisation  of 
F(f2).  To  achieve  this,  let  the  function  <j>  evolve  as  a  function  of  the  time  t  according  to  an  advection  equation 

=  -v(x,t)  ■  Vfax,t)  (2) 

where  v  is  a  local  velocity,  to  be  specified  below.  Assume  that  for  t  =  0  then  (j>  solves  the  eikonal  equation  |V0|  =  1. 
Given  that  (f>  =  0  on  the  boundary  T,  this  is  equivalent  to  taking  fax)  to  be  the  (signed)  distance  of  point  x  from 
this  boundary.  The  vector  V(f>  is  normal  to  the  level  sets  of  <f>,  and  it  is  clear  from  |2|  that  the  time-evolution  of  4> 
depends  only  on  the  normal  velocity  vn  =  v  ■  V<^>/|V0|.  In  what  follows  we  take  V«n  •  V^>  =  0  which  ensures  that  if 
the  eikonal  equation  is  true  at  t  =  0,  then  |V</>|  =  1  for  all  times  t  >  0. 

Hence,  to  specify  the  time  evolution  of  <j>,  it  is  sufficient  to  specify  vn  for  all  points  on  the  boundary  T  of  the  structure 
Q.  since  the  condition  Vu”  •  V0  =  0  then  specifies  u11  at  all  other  points.  In  order  that  ©  corresponds  to  steepest 
descent  for  the  objective  function  F,  we  introduce  a  sensitivity  sF ,  so  that  sF( u)  is  the  sensitivity  for  F  at  the  point 
X  (it).  To  define  sF ,  consider  a  deformed  structure  whose  boundary  Te  consists  of  points  Xe(u)  =  X (u)+ez(u)n(u) 
where  the  function  z  sets  the  size  of  the  displacement  of  each  boundary  point.  Informally,  the  sensitivity  sF (u)  is  the 
rate  of  change  of  F  associated  with  moving  the  boundary  point  X(u)  in  the  direction  n{u).  More  precisely,  sF  is  the 
unique  function  that  obeys 


Ffale)  =  F(ft)  +  e  j  sf{u)z(u)<M(u)  +  0(e2).  (3) 

We  assume  in  the  following  that  the  objective  function  F  and  the  boundary  T  are  sufficiently  smooth  that  this 
sensitivity  function  exists.  For  a  rigorous  discussion  of  these  issues,  see  [T] . 

To  perform  an  unconstrained  minimisation  of  W (SA) ,  one  should  prepare  an  initial  condition  in  which  </>  is  the  signed 
distance  from  the  boundary  of  some  initial  design  Ho-  Then  one  should  solve  ©  ,  with  the  normal  velocity  at  boundary 
point  X(u)  being 


«"(*(«))  =  sF(u).  (4) 

The  normal  velocity  vn  at  any  point  x  which  is  not  on  the  boundary  T  should  be  set  equal  to  the  normal  velocity  at 
the  nearest  point  on  T,  which  ensures  that  Vn11  •  V0  =  0,  as  noted  above.  This  time-evolution  for  <j)  encodes  a  time- 
evolution  for  the  structure  ft:  given  that  (j)  evolves  in  this  way,  it  is  easy  to  verify  that  dtF(Q)  =  —  f  sF  (u)2<M(u)  <  0, 
so  the  objective  function  decreases  with  time. 

We  now  consider  minimisation  of  F(fl)  subject  to  a  constraint,  that  G(fi)  >  0.  To  achive  this,  a  Lagrange  multiplier 
/rG  is  introduced  m,  such  that 


vn(X{u))  =  — sF  (n)  +  nGsG{u)  (5) 

where  sG  is  the  sensitivity  for  the  constraint  function  G ,  and  /uG  is  chosen  (independent  of  u )  such  that  the  system 
does  not  violate  the  constraint  on  G.  The  generalisation  of  this  construction  to  systems  with  multiple  constraints 
or  to  equality  constraints  of  the  form  H(Q)  =  0  is  straightforward  but  we  restrict  here  to  just  one  constraint,  for 
compactness  of  notation. 


B.  Discrete  space  and  time,  and  boundary  discretisation 

For  a  computational  implementation,  these  equations  must  be  discretised  in  both  space  and  time.  For  the  spatial 
discretisation,  the  domain  Sld  is  partitioned  into  a  square  grid  of  side  a0  =  1:  see  Fig.  [l]  The  vertices  of  the  grid 
are  called  nodes.  Each  node  i  has  a  position  x j  and  an  associated  value  of  the  level-set  function  fa  =  faxi).  The 
temporal  discretisation  is  a  simple  first-order  Euler  scheme,  so  ([2])  becomes 

fa{t  +  At)  =  fa{t)  -  vf(t)At  ■  |V0(i)|<  (6) 

where  v f  is  the  normal  velocity  at  node  i  and  \S7(f>(t)\i  is  the  modulus  of  V$,  evaluated  at  a which  is  equal  to  unity  if 
<f>  solves  the  eikonal  equation  exactly.  In  practice,  V0  is  estimated  for  each  node  using  the  Hamilton- Jacobi  weighted 
essentially  non-oscillatory  method  (HJ-WENO)  described  in  [IS]. 

In  order  to  determine  the  normal  velocities  {?;“},  the  LSTO  method  employs  a  second  level  of  discretisation:  see 
Fig.  i  From  the  (discretised)  level  set  function  <j>,  a  (discrete)  set  of  boundary  points  is  inferred,  as  follows:  If 
fa  =  0  then  node  i  is  a  boundary  point.  Also,  if  two  adjacent  nodes  have  ^-values  with  opposite  signs  then  there  is 
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FIG.  2.  Enlarged  view  of  the  discretised  zero  contour  of  the  level  set  function.  Each  boundary  point  a  is  associated 
with  two  segments  of  the  piecewise  linear  boundary.  For  each  boundary  point,  we  define  an  associated  boundary  length 

£a  =  {la,  1  +  h,a,2)  /2. 


a  boundary  point  between  them,  which  is  taken  to  lie  on  the  edge  between  the  nodes,  with  a  position  determined  by 
linear  interpolation.  Let  the  number  of  boundary  points  be  n  and  let  the  position  of  boundary  point  a  be  Xa  (with 
1  <  a  <  n).  The  boundary  points  form  a  set  of  closed  curves,  which  provide  a  discrete  representation  of  the  boundary 

r. 

The  reason  for  this  boundary  discretisation  is  that  the  LSTO  method  uses  estimates  of  the  sensitivities  sF,sG 
on  the  boundary  points  (this  is  the  natural  choice  since  the  sensitivity  is  intrinsically  related  to  the  boundary  T). 
Given  these  sensitivities,  one  infers  a  velocity  V™  for  each  boundary  point  a:  these  velocities  are  determined  by  a 
simple  optimisation  sub-problem  that  is  solved  at  each  time  step.  From  the  boundary  point  velocities  Vj",  the  normal 
velocities  v f  for  each  node  are  calculated  by  a  fast-marching  method.  Hence  the  level  set  variables  can  be  propagated 
forward  in  time,  according  to  ([6]).  We  now  describe  these  steps  in  more  detail. 


C.  Determination  of  boundary  point  velocities  from  an  optimisation  sub-problem 


To  determine  the  boundary  point  velocities  V£,  we  suppose  that  estimates  of  the  sensitivities  sF  and  sG  are 
available  for  each  boundary  point.  For  boundary  point  Xa  let  these  estimates  be  sF,sG.  (Estimation  of  these 
sensivities  depends  on  the  problem  of  interest  and  will  be  discussed  in  Sec.  |IV|)  Also  define  lai\  and  la ^  as  the 
distances  from  point  a  to  its  neighbouring  boundary  points;  and  let  £a  =  {la,i  +  la, 2)/^  be  the  length  of  boundary 
associated  with  point  a,  as  in  Fig.  [2]  (It  follows  that  is  the  total  boundary  length.)  Now  suppose  that  each 

boundary  point  Xa  moves  a  distance  V^At  in  the  normal  direction:  discretising  ([3])  along  the  boundary,  the  change 
in  the  objective  and  constraint  functions  can  be  estimated  as 

n 

=  E  Oa^Af, 

E  (7) 

AG  =  ^  V*sG£aAt. 

Oi= 1 


Direct  optimisation  of  the  V)"  can  then  be  used  to  optimise  the  change  in  F,  given  any  constraints  (see  for  example  [7]). 

However,  optimising  over  all  the  parameters  V™  is  not  convenient  numerically:  there  is  a  large  number  of  such 
parameters,  and  discretisation  errors  can  result  in  rough  boundaries  T  JT] .  Instead,  the  LSTO  method  uses  the 
(constrained)  steepest-descent  defined  by  J5| ,  with  a  variable  time  step  that  is  optimised  according  to  the  values  of 
the  sensitivities.  Spatial  discretisation  of  ([5^  yields 


WnA  t  =  X  FsF  +  Ac 


(8) 


where  \F ,  XG  are  parameters  to  be  determined  (with  —  XF  corresponding  to  the  time  step  At  in  ([5])  and  XG  corre¬ 
sponding  to  [iG At). 

At  each  iteration,  the  LSTO  method  optimises  the  parameters  XF,XG,  in  order  to  make  A F  as  negative  as  possible. 
However,  this  procedure  is  subject  to  several  constraints,  which  include  both  the  optimisation  constraint  ([I]),  and 
several  additional  considerations,  which  are  determined  by  the  physical  characteristics  of  the  topology  optmisation 
problem.  Note  that  the  handling  of  these  constraints  in  the  algorithm  presented  here  differs  from  the  method  of  [7]. 
The  method  presented  here  has  been  chosen  to  combine  simplicity  and  accuracy. 
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First,  note  that  A F,  AG  are  estimates  for  the  changes  in  F,  G,  which  are  accurate  only  if  the  boundary  point 
displacements  are  not  too  large.  To  avoid  very  large  values  of  the  A  parameters,  a  constraint  is  applied  to  each 
boundary  point  displacement: 


|V^Ai|  <  dcFL 


(9) 


Note  that  consistency  between  ([8])  and  ([5])  shows  that  —  XF  is  equal  to  the  timestep  At:  the  primary  role  of  this 
constraint  is  to  impose  a  sufficiently  small  time  step  for  the  Euler  discretisation  (|6|.  To  determine  an  appropriate 
value  for  dcFL,  we  use  the  Courant-Friedrichs-Levy  (CFL)  condition  associated  with  ([2]),  which  requires  that  dcFL  be 
smaller  than  the  grid  spacing.  Large  dcFL  leads  to  larger  time  steps  and  hence  faster  optimisation,  but  if  the  resulting 
convergence  histories  are  interpreted  as  as  approximate  solutions  to  ©  then  these  solutions  are  less  accurate  when 
dcFL  is  large.  The  specific  value  of  dcFL  depends  on  the  problem  of  interest. 

To  implement  the  CFL  constraint,  the  optimisation  domain  for  XF  is  |AF|  <  A^ax  with  XFax  =  —dCFL/( maxQ  \sF\). 
A  similar  restriction  is  applied  to  AG.  Since  these  constraints  are  applied  separately  to  \F  and  AG,  the  resulting 
solution  may  still  violate  the  constraint  ©:  if  this  happens  then  the  resulting  A  parameters  are  rescaled  by  a  factor 
rfcFL/(maxQ  |V^At|),  and  the  Vff  are  recalculated  using  ([8]),  so  that  ©  is  then  satisfied. 

Second,  if  the  design  at  time  t  is  fit,  the  inequality  constraint  Q  for  the  optimisation  requires  G(flt)  +  AG  >  G*. 
If  fit  does  not  satisfy  the  constraint  (G(flt)  <  G*)  then  the  optimiser  may  not  be  able  to  find  any  solution  for  which 
AG  is  large  enough  to  solve  the  constraint.  For  this  reason,  the  optimisation  is  performed  subject  to  a  modified 
constraint  AG  >  Go  where  Go  =  G*  —  G(flt)  if  a  solution  is  possible,  otherwise  a  smaller  value  for  Go  is  chosen.  (In 
practice,  the  maximum  possible  AG  is  calculated  by  considering  the  cases  XF'G  =  ±A(^:  if  the  maximal  possible 
value  AGmax  is  less  than  G*  —  G(flt)  then  it  is  likely  that  the  constraint  cannot  be  satisfied  so  we  take  Go  =  cAGmax, 
where  the  parameter  c  can  be  adusted  according  to  the  problem  of  interest.  Typically  we  take  c  ss  0.5.  Note  that 
AGmax  <  —G(flt)  usually  happens  only  during  the  initial  stages  of  optimisation,  so  the  way  that  this  case  is  handled 
does  not  typically  affect  the  convergence  of  the  method  to  its  final  solution.) 

Third,  note  that  boundary  points  should  not  move  outside  the  design  domain  fl^.  To  avoid  this  problem,  let  da  be 
the  signed  distance  of  boundary  point  a  from  the  boundary  of  fid-  (The  sign  of  da  is  positive  if  the  normal  vector  na 
points  towards  the  domain  boundary,  negative  otherwise.)  We  estimate  that  point  a  has  moved  outside  the  domain 
if  VffAt  >  da  >  0  or  VffAt  <  da  <  0  (this  is  an  approximation  since  the  normal  velocity  is  not  perpendicular  to  the 
interface,  but  it  is  sufficient  for  our  purposes).  If  this  happens  we  replace  ^  by  VffAt  =  da. 

With  these  ingredients  in  place,  we  finally  define  the  full  optimisation  problem  that  is  used  to  determine  XF , XG . 
Combining  ®,  and  accounting  for  the  possibility  that  boundary  points  might  move  outside  of  the  domain  Aid,  define 


AF(Af,Ag)  =  ^ 


'a  Za^ot , 


AG(Xf,Xg)  = 


(10) 


where  za  =  X FsF  +  A GsG  if  the  boundary  point  remains  inside  the  domain,  and  za  =  da  otherwise.  We  then  choose 
XF'G  to  minimise  AF(Xf,  Ag)  subject  to  AG(XF ,  AG)  >  Go,  on  the  domain  |AF|  <  A^)ax,  |AG|  <  AGax.  Given  the  XF’G 
solving  the  optimisation  problem,  the  time  step  is  At  =  —  XF  and  the  boundary  point  velocities  are  Vff  =  zaj At. 

In  practice,  this  two-parameter  optimisation  sub-problem  is  solved  using  the  SLSQP  method  from  the  NLOPT 
package  (251 .  The  simple  form  of  (10)  means  that  derivatives  of  AF ,  AG  can  be  obtained  analytically  without  the 
need  for  the  finite  differences  used  in  |7| .  If  some  sensitivities  are  very  large,  is  convenient  to  precondition  the  optimiser 
by  defining  rescaled  parameters  XF  =  A F  / A  and  corresponding  rescaled  sensitivities  sF  =  AsF  for  some  parameter 
A,  to  ensure  efficient  solution  of  this  sub-problem. 


D.  Level  set  update 

Having  calculated  the  boundary  point  velocities  Vff,  the  velocities  vf  on  nodes  adjacent  to  the  boundary  points  are 
calculated  by  inverse-squared  distance  weighting.  That  is,  if  the  edges  associated  with  node  i  of  the  grid  include  m 
boundary  points,  then  these  points  are  indexed  by  [3  =  1 ...  to.  Let  their  velocity  of  point  /3  be  Vff  and  its  distance 
from  node  i  be  rp  .  Then  vf  =  X!fci(^/rs)/E™=i(Vr/3)]-  This  fixes  the  uf  on  a  narrow  strip  that  contains  the 
boundary  T.  To  fix  the  vf  on  the  other  nodes  requires  a  velocity  extension  procedure,  for  which  we  use  a  fast-marching 
method,  as  in  [7]. 

Given  the  vf  on  all  nodes,  the  level  set  variables  </>j  are  updated  according  to  ([6]).  The  whole  process  -  inference 
of  boundary  points,  sensitivity  calculation,  optimisation  of  XFXG ,  velocity  extension,  level  set  update  -  is  repeated 
until  the  algorithm  converges  to  an  optimal  structure. 
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Finally,  note  that  in  practice,  it  is  convenient  to  do  the  velocity  extension  and  the  level  set  update  only  in  a 
narrow  band  close  to  the  boundary.  This  improves  the  efficiency  of  the  method  but  it  means  that  4>i  is  given  by 
the  signed  distance  to  the  boundary  only  within  the  narrow  band.  To  correct  for  this  effect,  all  of  the  4>i  variables 
are  periodically  reinitialised  to  be  consistent  with  a  signed  distance  function.  This  reinitialisation  uses  the  same 
fast-marching  implementation  used  for  the  velocity  extension. 

Note  that  given  a  set  of  boundary  points  Xa ,  reinitialisation  of  the  4>i  followed  by  a  recalculation  of  the  boundary 
point  positions  does  in  general  lead  to  small  changes  in  the  boundary  point  positions  (see  appendix).  As  the  system 
gets  close  to  convergence,  this  can  have  small  but  significant  effects  on  the  objective  function,  which  acts  as  a  weak 
source  of  numerical  noise  during  optimisation.  This  fact  has  implications  for  the  accuracy  of  the  stochastic  method 
described  below. 


III.  STOCHASTIC  LEVEL  SET 

Having  described  the  deterministic  optimisation  algorithm  of  izmsi,  the  next  step  is  to  introduce  a  stochastic 
component  to  this  algorithm.  Several  other  stochastic  level-set  methods  have  been  considered  recently  maun  mi, 
but  differ  from  the  approach  proposed  here  in  that  the  noise  in  this  scheme  is  applied  directly  only  on  the  boundary  T, 
which  ensures  that  the  function  <f>  retains  its  property  as  a  signed-distance  function  under  the  stochastic  time-evolution. 

We  first  consider  the  optimisation  problem  for  F(Sl)  in  the  absence  of  any  constraint  on  G(H).  In  this  case  the 
deterministic  algorithm  corresponds  to  steepest  descent  for  F,  which  converges  to  a  local  optimum.  By  contrast,  the 
corresponding  stochastic  algorithm  does  not  converge  to  an  optimum  of  F.  Instead  it  converges  to  a  steady  state 
in  which  it  explores  a  range  of  structures  H.  The  algorithm  is  designed  so  that  the  probability  pr  that  it  visits  a 
structure  fl  within  the  steady  state  is  proportional  to  a  Gibbs-Boltzmann  factor 

Pt(H)  oc  e-F(n^T  (11) 

where  T  is  a  noise  intensity  which  in  physics  would  be  identified  as  a  temperature.  For  T  — >  0  we  see  that  the 
probability  will  concentrate  close  to  the  optimal  structure  f2*,  which  minimises  F. 

The  advantage  of  the  stochastic  algorithm  is  that  it  can  explore  many  different  (non-optimal)  designs.  In  particular, 
for  a  non-convex  optimisation  problem  and  given  a  long  enough  time  t,  it  should  explore  all  the  local  optima,  not  just 
the  one  closest  to  the  starting  point  (which  would  be  found  by  a  deterministic  method).  Moreoever,  the  distribution 
of  structures  that  the  algorithm  finds  is  controlled  via  This  means  that  the  steady  states  of  the  algorithm  at 
different  temperatures  are  related,  as 


PT3(fl) 
PTi  (^) 


oc  e 


F(n)(i/T!-i/T2) 


and  in  particular,  the  marginal  distribution  of  the  objective  function  itself  satifies 

Pt2(F)  oc  eF(1/Tl-1/T2)PTl(F) 


(12) 


(13) 


Such  relationships  between  distributions  form  the  basis  of  simulated  annealing  and  parallel  tempering  methods  0111 
m ,  which  have  proven  useful  in  exploring  many  non-convex  optimisation  problems. 

Note  however  that  these  results  are  based  on  the  proportionality  relationship  More  generally,  we  expect  that 
the  invariant  measure  of  the  stochastic  process  for  the  structures  is 


dPT(tt)  =  z^T)e~F{n)/Td P°(n) 


(14) 


where  Z(T )  is  a  normalisation  constant  and  p°  is  a  reference  measure  for  structures  (independent  of  T).  We  do  not 
have  a  precise  characterisation  of  the  measure  p°,  but  as  long  as  p°  is  independent  of  T  then  ( 12|l3 1  are  valid,  and 
can  be  used  to  check  consistency  between  our  algorithm  and  the  asserted  invariant  measure  <E3b 


A.  Stochastic  motion  of  boundary  points 


To  describe  our  stochastic  method,  it  is  useful  to  cast  ©  as  an  equation  of  motion  for  the  boundary  points  Xa: 

d 


dt 


Xa(t)  =  V«(t)na(t). 


(15) 
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where  na  is  a  unit  vector  in  the  direction  of  the  (outward)  normal  to  the  boundary  T  at  point  Xa.  The  stochastic 
element  of  the  dynamics  operates  directly  on  the  boundary  points.  To  implement  this,  Eq.  ©  is  replaced  by  a 
stochastic  differential  equation  (SDE).  For  the  deterministic  problem  of  minimisation  of  F  in  the  absence  of  any 

U|.  Then  a  natural  generalisation  of  (15)  is  the  SDE: 


constraint,  we  write  V£  =  —  sF,  consistent  with 


dX“  =  -si 


,dt  +  nay/2T/£a  odWt° 


(16) 


The  theory  of  such  equations  is  discussed  (for  example)  in  ^ljjj.  Roughly  speaking,  one  may  interpret  dW“  as  a 
small  increment  in  the  boundary  point  position  Xa,  associated  with  a  small  time  increment  d t.  The  increment  dX“ 
consists  of  a  determinstic  part  —sFnadt  that  is  proportional  to  dt,  and  a  random  part  na  ^2T /£a  o  d The  o 
indicates  that  the  SDE  is  written  in  the  Stratonovich  convention:  the  implications  of  this  will  be  discussed  below.  The 
increment  dW“  is  a  standard  white  noise  (or  Wiener  process)  associated  with  boundary  point  a,  which  is  independent 
of  the  noises  on  all  other  boundary  points. 


Inspection  of  ( 16 )  shows  that  the  noise  term  acts  in  a  direction  perpendicular  to  the  boundary  (as  might  be 
expected);  it  is  proportional  to  y/T  which  sets  the  noise  strength.  The  factor  of  yjl/£a  might  not  be  expected  a  priori 
-  it  is  necessary  because  consistency  of  the  stochastic  level-set  method  requires  that  the  noise  intensity  is  equal  at  all 
points  on  the  boundary  T,  but  the  boundary  points  are  not  equally  spaced  along  T  (recall  Fig.  [2]) .  The  idea  is  that  a 
given  boundary  T  has  many  possible  discretisations  in  terms  of  boundary  points,  but  the  resulting  stochastic  evolution 
should  be  independent  of  this  discretisation.  Mathematically,  this  idea  can  be  encapsulated  as  a  reparameterisation 
invariance  of  the  equations  of  motion,  as  we  now  discuss. 


B.  Evolution  of  continuous  boundaries,  and  reparameterisation  invariance 


It  is  useful  to  consider  a  process  by  which  a  continuous  curve  evolves  in  time,  and  to  interpret  ( 16 )  as  a  discretised 
approximation  of  this  process.  For  the  continuous  curve  we  write  a  generalisation  of  Q,  as 


dXt(u)  =  —sF  (u)n(u)dt  +  n(u) yj2T / £' (u)  o  d Wt(u) 


(17) 


where  Wt(u)  a  random  Brownian  noise  with  equal  intensity  at  each  point  u.  (Technically,  the  quadratic  variation  of 


this  process  satisfies  d(/“2  Wt(u)du)  =  («2  —  Ui).)  To  derive  (16)  as  a  discretised  version  of  (17),  we  identify  Xa  as 


the  centre  of  mass  of  a  small  segment  of  the  boundary:  =  (l/£a)  Ju  2  Xt(u)d£(u).  Assuming  (see  below)  that  the 

segment  is  small  enough  that  n(u)  and  i'(u)  and  sF(u )  can  be  replaced  by  their  mean  values  within  the  integral,  this 
yields  (fl6|). 

Within  this  continuous  setting,  the  idea  that  the  evolution  of  the  curve  should  be  independent  of  its  discretisation  in 
terms  of  boundary  points  has  a  mathematical  statement  in  terms  of  reparameterisation  of  the  coordinate  u.  Consider 
a  closed  curve  described  by  described  by  the  internal  co-ordinate  u  £  [0, 1],  so  X  :  [0, 1]  — >  fid  is  a  smooth  function 
with  X(u)  a  point  on  the  curve.  Now  consider  a  continuous  monotonically-increasing  function  /  :  [0, 1]  — >  [0, 1]  with 
/( 0)  =  0  and  /( 1)  =  1.  Define  X  :  [0, 1]  — >  by  X(u)  =  X(f(u)).  It  should  be  clear  that  X  and  X  are  different 
representations  of  the  same  curve. 

In  order  that  the  evolution  of  a  curve  T  does  not  depend  on  its  parameterisation  in  terms  of  boundary  points,  we 
require  that  (17)  evolves  the  functions  X  and  X  in  the  same  way.  This  may  be  verified  by  direct  substitution,  as 
long  as  the  noise  prefactor  y/ljtjuj  is  included.  This  is  the  reason  for  including  the  factor  in  y/l/£a  that  multiplies 


the  noise  in  ( 16 ). 


Note  that  we  have  assumed  throughout  that  the  boundary  T  is  smooth  enough  that  the  normal  vector  n(u)  can  be 
defined  at  every  point  X  (u) .  Even  for  deterministic  optimisation,  this  assumption  requires  some  care,  if  the  optimal 
shape  has  a  boundary  that  includes  kinks.  For  stochastic  optimisation,  there  is  an  additional  factor,  which  is  that 
the  noise  tends  to  roughen  the  curve  and  the  normal  vector  n(u )  may  not  be  defined.  In  our  numerical  scheme, 
the  boundary  is  discretised  and  a  WENO  estimate  of  V<(>  can  always  be  used  to  define  a  local  normal,  so  potential 
problems  with  rough  boundaries  do  not  appear.  For  a  detailed  mathematical  analysis,  some  regularisation  is  required 
to  ensure  existence  and  uniqueness  of  the  SDE  ©,  but  we  defer  these  issues  to  a  later  publication. 


C.  Invariant  measure 


Our  assertion  above  was  that  the  invariant  measure  for  this  process  is  proportional  to  (Tl]).  We  do  not  have  a 
rigorous  proof  of  this  claim,  although  we  will  demonstrate  below  that  our  numerical  method  is  consistent  with  (13). 
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FIG.  3.  Shape  matching  by  the  level-set  method  (without  any  stochastic  component).  The  target  shape  is  a  discretised 
representation  of  the  Stanford  bunny  (top  right);  the  initial  shape  is  a  circle  (top  left).  The  main  figure  shows  a  time  series 
of  the  (normalised)  objective  function  .F(fit)/A(fid)  with  F  given  by  (|24|)  and  A(fid)  the  total  area  of  the  design  domain. 
Representative  shapes  obtained  during  this  convergence  proedure  are  also  shown  (above).  The  inset  illustrates  the  definition 
of  the  objective  and  sensitivity  functions,  for  the  simpler  case  of  a  circular  target  shape  and  and  a  square  initial  shape  fi. 
Blue  shading  indicates  regions  where  <)>*arset  <  fa  and  sF  =  —  1,  so  outward  motion  of  the  interface  reduces  F.  Similarly  green 
regions  have  </rarget  >  fa  and  sF  =  +1.  The  normalised  objective  function  .F(fi)/A(fid)  is  given  by  the  sum  of  the  blue  and 
green  areas,  as  a  fraction  of  the  total  area  of  the  design  domain  fid- 


As  a  plausibility  argument  for  our  scheme,  we  note  that  if  nasF  in  (16)  is  equal  to  (dF/3XQ)£a  1  then  (16)  has  the 
form 


dX“  =  ( 8F/dXa)a2dt  +  aV2T  o  dW° 


(18) 


with  a  =  1  / v^q ■  It  is  well-known  [14]  that  the  invariant  measure  for  this  equation  is  of  the  form  (11).  Moreover, 
spatially  discretising  ([3])  and  identifying  ez(u)  as  the  boundary  point  displacement  na  ■  5Xa  yields 


F( fie)  -  F(Q)  =  ^2  sFiana  ■  bXa  +  0(SX2) 


(19) 


indicating  that  indeed  £anasF  does  indeed  correspond  to  the  derivative  (dF/dXa),  as  required.  Refining  this  argu¬ 
ment  into  a  more  rigorous  proof  is  an  interesting  direction  for  future  work. 

In  the  following,  we  show  how  the  process  ( 16 )  can  be  implemented  within  the  level-set  method  described  in  Sec.  [II] 
We  assume  that  ( 12 1  holds  theoretically  and  we  investigate  the  extent  to  which  it  holds  numerically. 


D.  Stochastic  dynamics  with  a  finite  time  step 


To  implement  the  stochastic  evolution  (16)  within  the  level-set  method  of  Sec.[TTJ  we  require  a  generalisation  of  @. 
To  achieve  this,  integrate  ( 16 )  over  a  small  time  interval  At  and  identify  the  average  boundary  point  velocity  in  the 
normal  direction  as  V£(t)  =  (At)-1  ra(s)  •  dX“.  Hence,  keeping  terms  up  to  first-order  in  At: 


V£(t)At  =  — *£(t)At  + 


(20) 


where  £  is  a  standard  normal-distributed  random  number.  The  first  term  on  the  right  hand  side  is  the  standard 
deterministic  increment  corresponding  to  (15 1.  The  second  term  is  a  random  increment  (with  zero  mean).  Note 
however  that  the  length  ta  in  this  second  term  is  evaluated  at  time  t  +  At/2,  which  corresponds  to  the  midpoint  of 
the  time  interval  (this  is  the  meaning  of  the  Stratonovicli  product  that  was  denoted  by  o  in  @).  To  arrive  at  an 
equation  that  can  be  used  in  a  (first-order)  numerical  scheme,  it  is  necessary  to  calculate  V£(t)  in  terms  of  quantities 
that  are  evaluated  only  at  time  t.  This  can  be  achieved  by  using  Ito’s  formula  [14],  which  yields 


Vs{i)At  =  -sFa{t)At  + 


1 2TAt  t  Tk0 

2 Ta 


At. 


(21) 
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(c) 


P(F) 


(d) 


FIG.  4.  Shape  matching  with  finite  noise  strength  T:  (a)  Time  series  of  the  area  mismatch  (24 1  for  several  noise  strengths 
(temperatures).  Raising  temperature  leads  to  increasingly  sub-optimal  designs  and  larger  fluctuations,  (b)  A  representative 
shape  (red)  from  the  steady  state  at  temperature  T  =  1,  compared  with  the  target  shape  (red).  Mismatch  with  the  target 
structure  is  most  pronounced  in  regions  of  high  curvature,  (c)  Histograms  (normalised  as  probability  densities)  of  the  objective 
function,  obtained  from  long  time  series  in  the  steady  state  of  the  system.  Data  are  shown  for  seven  temperatures  equally 
spaced  between  T  =  0.1  and  T  =  0.25,  with  the  remaining  temperatures  equally  spaced  between  T  =  0.25  and  T  =  1.0.  (d)  The 
measured  probability  densities  agree  well  with  the  predictions  of  Eq.  ( 13 1 :  we  show  measured  pr{F)  for  various  temperatures, 
and  the  corresponding  predictions  using  Eq.  (|13|)  with  Ti  =  T  and  T2  =  T  +  AT.  The  temperature  increment  AT  is  such  that 
the  distributions  at  Ti,T 2  are  adjacent  in  (c). 


where  na  is  the  (signed)  curvature  of  the  boundary  at  point  a  (see  Appendix). 

The  derivation  of  ©  from  (p~6])  is  straightforward  within  the  framework  of  stochastic  calculus.  We  omit  technical 
details  and  provide  a  short  argument  to  justify  it:  for  a  boundary  point  increment  dXa,  the  change  in  the  length 
of  the  boundary  segment  is  dla  =  na  •  («;„  o  dXQ).  (Roughly  speaking  this  corresponds  to  the  chain  rule,  with 
d£a/dXa  =  naKa,  as  discussed  in  the  Appendix.)  Ito’s  formula  [14]  states  that  for  an  SDE  of  the  form  dxt  = 
f(xt)dt  +  cr(xt)  °  d Wt,  the  appropriate  first-order  discretisation  is  Ax  =  f(xt)At  +  a(xt)^V~Ai  +  \a' (xt)a(xt)At.  In 
this  case  we  have  from  (16)  that  a  =  na\/2T/£a:  the  analogue  of  a'  requires  a  derivative  with  respect  to  £ai  which 
yields  (—nan/2la)^/2T/£ui  and  hence  (21). 


E.  Implementation  in  level-set  method,  and  incorporation  of  constraints 


The  stochastic  velocity  (21)  is  straightforwardly  incorporated  into  the  deterministic  level  set  method  described  in 
Sec.  [TT]  This  yields  our  stochastic  level  set  optimisation  method.  The  only  modification  of  the  deterministic  algorithm 
is  that  noise-dependent  terms  are  added  to  We  consider  the  case  of  optimisation  subject  to  a  constraint  G  =  0 
as  in  the  deterministic  case. 

In  the  stochastic  method,  the  parameters  XF ,  XG  in  |8)  are  calculated  exactly  as  in  the  deterministic  method.  The 
simplest  approach  is  then  to  identify  XF  — >  — At  and  add  the  stochastic  terms  in  (21)  to  the  boundary  increments. 
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FIG.  5.  Time  series  for  deterministic  minimisation  of  shape  perimeter  (circles,  left  axis),  subject  to  a  constraint  on  that  the 
area  fraction  outside  Q  must  be  at  most  60%  of  the  total  area  (triangles,  right  axis).  The  top  panel  shows  the  time-dependence 
of  Q,  starting  from  a  square  initial  shape,  which  violates  the  constraint.  The  shape  S2  increases  in  size  (to  satisfy  the  constraint) 
and  evolves  towards  a  circle  (to  minimise  the  perimeter).  The  domain  fid  is  a  square  grid  of  size  200  x  200. 


However,  since  the  first  stochastic  term  is  proportional  to  V At,  this  can  lead  to  large  boundary  point  displacements 
that  violate  the  CFL  condition.  Hence,  given  a  deterministic  time  step  At  =  — AF,  we  calculate  a  typical  stochastic 
increment  Acctyp  =  V%TAt.  If  Aa;typ  >  dcFL/2,  the  deterministic  parameters  XF,XG  are  rescaled  by  a  factor 
dcFL/2Aa;typ.  That  is 


^stoch  ^  nrin(l,  dcFL/2A:rtyp)  (22) 

^stoch  ^  min(l,  dcFL/2AcctyP). 

This  construction  amounts  to  fixing  a  maximal  time  step  that  ensures  that  the  CFL  condition  is  obeyed  even  in  the 
presence  of  the  noise,  as  long  as  none  of  the  £a  parameters  are  too  small  (compared  to  unity).  (We  also  tested  a 
method  where  the  Aa;typ  is  adjusted  to  account  for  the  possibilities  that  some  £as  are  very  small,  but  this  led  to  a  less 
efficient  method  and  had  very  little  effect  on  the  results.)  The  time  step  is  — A^och  and  the  boundary  point  increment 
is  finally 


V£At  =  s 


FXF 

a  /'stoch 


+ 


Tna 

2£a 


X 


F 

stoch 


,  \G  G 

'  /xstochc’a! 


(23) 


This  equation  replaces  ^  within  our  stochastic  level  set  method.  Once  the  boundary  point  velocities  V£  have  been 
calculated  in  this  way,  the  rest  of  the  method  follows  exactly  the  deterministic  case:  the  velocities  on  the  nodes  are 
calculated  by  interpolation  and  fast-marching,  and  the  level  set  variables  are  updated  using  ([6]). 


IV.  RESULTS 

In  the  following  we  consider  several  examples  of  the  stochastic  level  set  method,  with  increasing  levels  of  complexity 


A.  Shape  matching 

We  first  consider  a  simple  geometric  optimisation  problem.  The  idea,  illustrated  in  Fig.  [3]  is  that  the  domain  Q 
should  match  a  predefined  reference  shape  retarget ,  in  this  case  a  discretised  two-dimensional  representation  of  the 
Stanford  bunny.  The  objective  function  is 


*■(«)  =  £  ^target  _j4.(Q) | 
i 


(24) 
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FIG.  6.  Stochastic  optimisation  of  perimeter  with  a  constraint  on  the  enclosed  area.  (a,b)  Time  series  for  the  objective  function 
F  (perimeter)  and  the  constraint  function  G  (area).  The  legend  in  (a)  applies  to  both  panels,  (c)  The  zero-temperature  optimum 
(a  circle),  compared  with  a  representative  sample  from  a  stochastic  calculation  at  T  =  0.5.  (d)  Histograms  of  the  objective 
function,  taken  from  the  steady  state  of  the  stochastic  process,  at  various  temperatures,  (e)  Measured  probability  densities, 
compared  with  the  predictions  of  Eq.  (13 1.  The  predictions  use  Eq.  (13 1  with  Ti  =  T  and  T2  =  T  +  AT,  as  in  Fig.  [4]  The 
agreement  is  good. 


where  the  sum  runs  over  the  cells  of  the  grid,  A,;(f2)  is  the  overlap  area  between  the  grid  cell  i  and  Q,  and  similarly 
^target  overiap  area  between  element  i  and  ^target- 

To  calculate  sensitivities,  note  that  if  A  is  the  area  enclosed  by  then  =  1  for  all  boundary  points  a.  For 
the  objective  function  (24),  this  requires  =  1  if  boundary  point  a  is  inside  retarget  and  si  =  —  1  otherwise.  To 
estimate  this  sensitivity  based  on  local  information,  we  calculate  for  each  node  i  its  squared  distrance  <^argelj  from  the 
boundary  of  the  target  shape.  For  each  boundary  point  a ,  we  estimate  its  signed  distance  from  the  boundary  of  the 
target  as  4>^arget,  which  is  interpolated  from  the  </ualgot  values  on  the  four  nearest  nodes.  We  also  interpolate  a  value 
(jjtrmi  ky  appjyirig  the  same  method  to  the  nearest  values  of  </>,:  since  the  boundary  point  is  on  the  zero  contour  of  (f> 
we  expect  |4Wlal|  <C  1  for  all  a.  The  sensitivity  is  then  estimated  as 


si  =  sign($taarget  -  $“al)  (25) 

Since  <jWlal  is  small  in  magnitude,  its  inclusion  in  |25|)  has  a  small  effect  when  the  current  and  target  shapes  are 
different  from  each  other,  but  ensures  that  the  deterministic  method  converges  (exactly)  to  12  =  £2tarset. 

Based  on  these  sensitivities,  Fig.[3]shows  the  deterministic  optimisation  algorithm  in  operation.  An  initially  circular 
structure  evolves  in  time  until  it  matches  exactly  the  target  shape.  The  CFL  constraint  is  dcFL  =  0.1  and  the  domain 
is  partitioned  into  a  grid  with  200  x  200  elements,  with  each  element  having  size  1.  The  initial  shape  is  a  circle  of 
radius  50. 

Fig.  a  shows  results  for  the  stochastic  algorithm,  applied  to  the  same  problem.  Panel  (a)  shows  that  the  objective 
function  does  not  converge  to  its  optimal  value  (F  =  0).  Instead  the  system  converges  to  a  steady  state  in  which 
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F  fluctuates  around  a  non-zero  mean  value.  As  the  noise  strength  increases,  both  the  mean  value  of  F  and  its 
fluctuations  increase,  as  the  algorithm  explores  a  range  of  non-optimal  designs.  Panel  (b)  shows  a  representative 
non-optimal  shape.  Panel  (c)  shows  distributions  (histograms)  of  the  objective  function,  for  different  values  of  the 
noise.  Finally  in  panel  (d)  we  test  the  prediction  of  (13),  that  distributions  of  the  objective  function  at  different  noise 
strengths  can  be  related  to  each  other.  That  is,  for  various  temperatures  T,  the  distribution  P(F)  can  be  predicted 
based  on  its  distribution  at  a  different  temperature  T  +  AT.  (This  prediction  is  possible  only  for  small  AT  since 
otherwise  the  exponential  factor  in  (13)  leads  to  large  statistical  errors.) 

The  results  shown  in  Fig.  Bd)  represent  strong  evidence  that  the  system  has  converged  to  a  steady  state  in  which 
the  predictions  of  ( 13 )  are  valid.  The  stochastic  algorithm  was  designed  in  order  to  obtain  a  steady  state  that  satisfies 
©,  and  the  numerical  agreement  with  (|13[)  is  a  fairly  stringent  test  of  this  criterion,  which  indicates  that  the 
theoretical  analysis  and  numerical  implementation  presented  here  are  self-consistent. 

Moreoever,  the  nature  of  the  non-optimal  shapes  shown  in  Fig.[4[b)  reveal  some  information  about  the  null  measure 
Po  in  (14):  of  the  shapes  12  for  which  F(fl)  is  not  equal  to  its  optimal  value,  the  method  seems  to  sample  preferentially 
those  shapes  with  low  curvature.  In  this  sense  the  numerical  noise  seems  to  act  to  smooth  the  boundary  T.  A  detailed 
analysis  of  this  effect  will  require  a  deeper  understanding  of  the  stochastic  dynamics  of  the  boundary,  which  is  an 
interesting  direction  for  future  work. 


B.  Perimeter  minimisation  with  area  constraint 


As  a  simple  optimisation  problem  that  includes  a  constraint,  we  next  consider  minimisation  of  the  perimeter  of 
f 2,  subject  to  a  constraint  on  its  total  area.  Given  that  F(Cl)  is  the  perimeter  of  f 2,  the  sensitivity  sF  is  is  given  by 
the  (signed)  Euclidean  curvature  n  of  the  boundary  T  (see  Appendix).  The  constraint  function  G  is  the  difference 
between  the  total  domain  size  and  area  of  f2:  that  is,  G(I2)  =  A(f2<j)  —  A(f2).  (We  consider  12  as  a  hole  in  large  piece 
of  material,  so  G(I2)  is  the  area  of  the  material.)  The  associated  sensitivity  is  sG  =  —  1. 

Estimation  of  the  curvature  k  is  non-trivial  for  the  discretised  boundaries  considered  here  [3].  Following  numerical 
tests  of  several  different  methods,  we  use  a  scheme  that  combines  accuracy  and  simplicity,  by  calculating  curvatures 
using  an  explicit  finite-difference  sensitivity  calculation.  Details  of  this  procedure,  and  its  associated  errors  are 
discussed  in  the  Appendix. 

Results  for  deterministic  optimisation  are  shown  in  Fig.  [5]  The  grid  is  of  size  200  x  200.  An  initially  square  structure 
12  evolves  into  a  circle  whose  area  satisfies  the  constraint.  In  this  case  G  >  0.6A(12d),  so  the  material  covers  at  most 
60%  of  the  domain  12^,  and  the  area  of  12  must  be  at  least  40%  of  this  domain. 

Fig.[6]shows  results  for  the  stochastic  method,  which  are  comparable  with  Fig. El  Fig.[6[b)  shows  that  the  stochastic 
algorithm  still  satisfies  the  constraint  (up  to  some  numerical  uncertainties).  Figf|6[e)  shows  that  the  results  are  again 
consistent  with  the  algorithm  sampling  a  distribution  of  shapes  consistent  with  (jnpl)- 

However,  the  example  shape  shown  in  Fig.  |6]c)  for  T  =  0.5  reveals  that  this  numerical  method  is  affected  by 
discretisation  effects  from  the  underlying  grid.  In  particular,  the  shape  12  shown  in  that  figure  is  not  circular,  but 
is  elongated  along  the  lattice  axes,  forming  a  kind  of  diamond  shape.  We  find  that  the  shapes  12  found  for  T  >  0 
consistently  have  this  property  (data  not  shown).  We  have  investigated  the  reason  for  this  effect,  which  we  attribute  to 
uncertainties  in  our  numerical  estimates  of  the  sensitivity  parameters  sa,  particularly  the  sensitivity  of  the  perimeter, 
which  is  the  curvature.  These  numerical  issues  are  discussed  in  the  Appendix.  In  terms  of  the  general  method,  our 
conclusion  is  that  the  stochastic  level-set  method  relies  on  accurate  sensitivity  estimates,  which  may  be  difficult  to 
obtain  in  practice.  However,  we  believe  that  the  method  itself  is  valid. 


C.  Non-convex  shape  matching  problem  with  perimeter  constraint 

So  far,  we  have  focussed  on  very  simple  optimisation  problems,  as  proof  of  principle  for  the  method.  However,  a 
central  motivation  for  the  development  of  a  stochastic  method  is  that  the  noise  forces  can  allow  the  system  to  visit 
multiple  locally-optimal  designs  in  a  non-convex  optimisation  problem.  As  a  very  simple  example  of  such  a  problem, 
we  consider  the  problem  shown  in  Fig.  [7] 

The  constraint  function  in  this  case  is  is  related  to  the  overlap  between  the  shape  12  and  the  dumbbell  (or  hourglass) 
target  shape  shown  in  Fig.  [7]  By  analogy  with  ([24]),  we  define 

G  =  ^|A‘arget-A4(12)|  (26) 

i 

where  the  sum  runs  over  all  cells  of  the  underlying  grid.  The  constraint  is  that  G  <  0.2A(12d):  that  is,  the  mismatch 
between  the  shape  12  and  the  target  (dumbbell)  shape  can  be  at  most  20%  of  the  total  design  domain.  In  practical 
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FIG.  7.  Non-convex  optimisation  problem.  (Top  row)  snapshots  from  an  optimisation  trajectory  in  which  the  perimeter  of 
the  red  shape  12  and  its  height  within  the  design  domain  are  both  being  minimised,  subject  to  a  constraint  on  the  mismatch  of 
the  red  shape  with  the  blue  “dumbbell’.  The  global  optimum  is  shown  in  the  rightmost  panel  while  the  second  panel  illustrates 
an  additional  local  optimum.  (Main  panel)  Time  series  of  the  objective  function  for  deterministic  optimisation  (T  =  0)  and 
stochastic  optimisation  (T  —  0.002).  Only  in  the  stochastic  case  does  the  system  escape  from  the  local  optimum  and  converge 
to  the  global  one.  The  mesh  size  is  100  x  100,  the  dumbell  consists  of  two  circles  of  radius  20  whose  centres  are  separated  by 
B  =  38.  The  parameter  a  =  0.65. 


terms,  this  means  that  when  f2  is  contained  within  one  of  the  lobes  of  the  dumbbell,  the  area  of  the  dumbbell  that  is 
outside  f2  must  be  at  most  0.2A(f2d).  The  dumbbell  is  given  by  the  union  of  two  circles  of  radius  20.  In  Cartesian 
coordinates,  the  size  of  the  design  domain  is  Lx  x  Ly  and  the  centres  of  the  circles  are  at  (Lx/ 2,  ( Ly  ±  B)j 2)  with 
B  =  38. 

The  objective  function  in  this  case  is  a  weighted  perimeter  F(f2)  =  M(fi)  with 


M(Sl)  =  /  m(X(u))<M(u) 


(27) 


If  m(X)  =  1  for  all  X  then  M  is  simply  the  perimeter  of  12.  The  idea  is  that  m(X)  is  the  weight  (per  unit  length) 
of  a  boundary  segment  at  position  X.  We  write  X  =  ( X ,  Y)  in  Cartesian  coordinates  and  rri(X)  depends  only  on  Y . 
It  is  given  by  to  =  1  for  positions  above  the  centre  of  the  upper  lobe  of  the  dumbbell,  and  m  =  a  for  positions  below 
the  centre  of  the  lower  lobe.  In  between,  m  varies  linearly  with  Y .  so  that 


i(X,Y)  = 


1, 

a, 

1  +  Q!  |  l  —  Oi 


l±z  +  ±F(Y-Ly/  2), 


Y  >  {Ly+B)/ 2 

Y  <  (Ly  -  B)/2 
otherwise. 


(28) 


The  parameter  a  =  0.65.  Physically,  the  objective  function  is  small  when  the  boundary  T  is  located  in  the  lower  lobe 
of  the  dumbbell,  and  larger  when  it  is  in  the  upper  lobe.  The  optimal  design  is  a  circle  located  inside  the  dumbbell, 

below  the  centre  of  the  lower  lobe.  _ 

The  sensitivity  sG  for  the  constraint  function  G  was  described  already  in  Sec.  IV  A  The  sensitivity  for  M  is 
=  Kam,(Xa)  +  nva  drnQYa ^  i  where  Ka  is  the  signed  curvature  and  nva  is  the  y  component  of  the  normal  vector  na. 
Fig.[?]shows  results  for  both  deterministic  and  stochastic  optimisation  for  this  problem.  The  determinstic  algorithm 
reveals  that  this  is  indeed  a  non-convex  optimisation  problem:  that  is,  there  are  two  local  minima  of  the  objective 
function.  Starting  with  a  circular  shape  located  in  the  upper  lobe  of  the  dumbbell,  the  deterministic  algorithm 
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FIG.  8.  Stochastic  optimisation  of  compliance.  (a,b)  Time  series  for  the  objective  function  (strain  energy)  and  constraint  func¬ 
tion  (material  area)  for  the  model  problem  described  in  the  text,  (c)  The  (local)  minimum  found  by  deterministic  optimisation 
(T  =  0)  and  a  representative  structure  at  T  =  0.0003.  (d)  Histograms  for  the  objective  function  for  different  noise  strengths 
T.  (e)  Test  of  the  reweighting  formula  (|13|).  The  results  are  not  consistent  with  (|13|) :  we  show  measured  distributions  at 


Tj  =  0.00020  and  Ti  =  0.00022.  The  distribution  at  T2  is  used  together  in  ( 13 1  to  arrive  at  a  prediction  for  the  distribution  at 
T\ ,  but  this  prediction  is  not  accurate.  We  attribute  this  effect  to  numerical  errors  associated  with  discretisation  and  sensitivity 
estimation,  as  discussed  in  the  main  text 


converges  quickly  to  a  local  optimum  in  which  the  shape  is  located  close  to  the  neck  (2nd  image  in  top  row  of  Fig.  [7|. 
The  global  optimum  has  the  shape  inside  the  lower  lobe,  but  convergence  to  this  shape  is  frustrated  because  the 
near-circular  locally-optimal  shape  cannot  fit  through  the  neck  of  the  dumbbell.  In  order  to  pass  through  the  neck, 
while  still  obeying  the  constraint,  the  area  inside  the  dumbbell  must  expand,  to  compensate  the  extra  area  outside. 
However,  this  leads  to  an  increase  in  the  objective  function  F,  which  is  not  possible  within  the  deterministic  algorithm. 

However,  on  adding  a  weak  stochastic  element,  one  sees  (for  these  parameters)  that  the  system  escapes  the  local 
optimum  and  converges  to  a  steady  state  where  it  samples  shapes  that  are  close  to  the  global  optimum.  A  final 
deterministic  optimisation  could  be  used  to  locate  the  true  optimum,  if  required.  This  example  shows  the  potential 
usefulness  of  the  stochastic  level-set  optimisation  method,  although  this  is  obviously  a  very  simple  model  at  this  stage. 


D.  Compliance  minimisation 

As  a  final  example,  we  show  how  this  method  can  be  applied  to  problems  inspired  by  engineering  applications  of 
shape  optimisation  (this  is  distinct  from  topology  optimisiation  in  that  no  holes  are  created  during  optimisation).  As 
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FIG.  9.  Approximations  for  the  local  Euclidean  curvature  around  a  boundary  point,  (a)  A  simple  geometric  approximation 
can  be  obtained  by  matching  a  circle  to  any  three  sequential  boundary  points.  The  three  points  are  assumed  to  lie  on  the 
perimeter  of  the  circle,  with  the  curvature  defined  as  the  inverse  of  the  circle’s  radius  [4].  (b)  The  curvature  can  also  be 
calculated  by  performing  an  explicity  finite-difference  sensitivity  calculation.  A  boundary  point  is  displaced  in  the  outward  and 
inwards  directions  along  its  normal  vector  and  a  sensitivity  is  defined  as  the  rate  of  change  of  the  discretised  perimeter  per  unit 
length,  i.e.  using  a  central  finite-difference.  In  the  limit  of  5  — >  0  this  defines  a  local  perimeter  sensitivity  for  the  boundary 
point,  (c)  Mean  Euclidean  curvature  for  circles  of  increasing  radius.  Red  circles  show  the  analytical  curvature,  1  / R.  Curvature 
measured  using  the  simple  geometric  approximation  (green  squares)  and  an  explicit  finite-difference  sensitivity  calculation  (blue 
triangles)  show  excellent  agreement.  A  value  of  5  =  IIP4  was  used  for  the  central  difference  calculation. 


a  proof  of  principle,  we  study  a  well-known  benchmark  problem  of  a  two-dimensional  cantilevered  beam  S3],  as  shown 
in  Fig.  [8]  The  overall  strain  energy  of  the  structure  is  minimised  subject  to  a  constraint  on  its  area.  The  sensitivity 
for  the  strain  energy  is  computed  using  a  linear  elastic  finite  element  analysis  and  the  sensitivity  associated  with  the 
area  constraint  is  sG  =  1,  as  above.  Once  the  sensitivities  are  known,  the  method  proceeds  exactly  as  in  the  simple 
examples  considered  in  previous  sections. 

In  this  problem,  the  grid  used  within  the  level-set  method  plays  a  second  role  as  a  finite-element  mesh  for  the 
structure.  The  objective  function  F  is  the  strain  energy  of  the  structure.  To  calculate  this,  the  structure  is  clamped 
at  the  left  boundary  ( Lx  =  0)  and  a  unit  load  is  applied  in  the  y-direction  at  position  (Lx,Ly/ 2).  The  Young’s 
modulus  of  the  material  inside  Cl  is  Y  =  100,  the  Poisson’s  ratio  is  0.3.  The  space  outside  Cl  is  occupied  by  a  weak 
material  whose  density  is  10~3.  The  strain  energy  is  minimised  subject  to  the  constraint  that  the  total  area  of  Cl 
is  at  least  0.5  of  the  domain  Cld-  The  sensitivity  for  the  compliance  is  calculated  by  a  standard  method  jTJ  and  the 
sensitivity  for  the  constraint  is  sG  =  +1,  as  described  above.  The  mesh  size  is  Lx  x  Ly  with  ( Lx,Ly )  =  (40,20). 
This  relatively  coarse  mesh  is  used  simply  for  computational  convenience  future  work  will  exploit  more  efficient 
sensitivity  calculations  which  will  allow  access  to  finer  grids,  but  this  is  beyond  the  scope  of  this  work. 

The  results  in  Fig.  [8]  are  qualitatively  consistent  with  Fig.  [6]  and  show  that  the  stochastic  level  set  method  can 
be  applied  in  such  contexts.  The  initial  condition  for  the  optimisation  is  a  completely  full  grid  Cl  =  Cld .  One  sees 
from  Figs.  |8](a,b)  that  the  optimiser  first  reduces  the  area  of  Cl  in  order  to  satisfy  the  constraint:  during  this  part 
of  the  algorithm  the  strain  energy  increases.  Then,  once  the  structure  satisfies  the  constraint,  its  shape  is  optimised 
to  reduce  the  strain  energy.  Fig.  [8](c)  shows  the  optimal  design  that  was  found  by  deterministic  optimisation,  and  a 
design  found  by  the  stochastic  method.  Fig.  [8](d)  shows  that  increasing  the  noise  increases  the  typical  values  of  the 
objective  function  F,  and  the  variance  of  this  quantity  also  increases.  However,  Fig.  [8j(e)  shows  that  the  results  are 
not  quantitatively  consistent  with  There  are  several  potential  sources  of  numerical  error  in  this  algorithm,  of 

which  the  largest  is  expected  to  be  the  numerical  uncertainties  in  calculations  of  sensitivities,  due  to  the  finite-element 
approximation.  We  attribute  the  deviations  from  the  predictions  of  ( 13 )  to  these  numerical  uncertainties. 

To  understand  this  effect  in  more  detail,  it  is  useful  to  notice  that  (13)  implies  quite  generally  that  Var(F’)  = 
T2^p(F),  where  ( F )  is  the  mean  strain  energy  at  temperature  T,  and  Var(F)  =  ( F 2)  —  ( F )2  is  its  variance.  This  is 
an  example  of  a  fluctuation-dissipation  theorem  (FDT)  [5].  For  the  systems  considered  in  Figs.  4j6  we  find  that  this 
relation  holds  (otherwise  the  predictions  based  on  the  reweighting  formula  would  fail).  For  this  strain  energy 
problem,  we  find  that  Var(i?)  «  3 T2-^p(F),  indicating  that  the  variance  of  the  objective  is  around  three  times  larger 
than  the  prediction  given  by  the  FDT.  Our  hypothesis  is  that  these  extra  fluctuations  in  F  come  from  numerical 
errors  associated  with  discretisation  and  with  estimates  of  the  sensitivity,  but  this  remains  an  area  for  future  study. 

Finally  we  note  that  the  globally-optimal  structures  for  this  (discretised)  compliance  problem  are  not  simply- 
connected  like  the  shapes  considered  in  Fig.  |8jc) :  the  optimal  structures  include  holes  [7j .  The  stochastic  method 
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FIG.  10.  Discretisation  errors  affect  mean  and  local  boundary  point  curvature,  (a)  Signed  distance  reinitialisation  using 
a  second-order  Fast-Marching  Method  introduces  errors  into  the  mean  boundary  point  curvature  (blue  squares).  When  the 
nodes  of  the  level  set  domain  are  initialised  using  an  exact  signed-distance  function  (using  the  exact  Euclidean  distance  to 
the  interface)  the  mean  curvature  agrees  near  perfectly  with  the  analytical  result  (red  circles).  In  this  case  the  only  errors 
are  introduced  by  the  piece-wise  linear  discretisation  of  the  boundary,  (b)  Local  boundary  point  curvature  around  a  circle 
of  radius  R  =  40.  The  angle  is  measured  relative  to  the  top  of  the  circle.  While  the  mean  curvature  is  excellent,  noise  is 
present  in  the  local  boundary  curvature,  even  when  using  a  perfect  signed-distance  function  (red  circles).  Reinitialisation  of 
the  signed-distance  function  leads  to  a  significant  increase  in  the  noise  (blue  squares). 


described  here  does  not  include  an  explicit  prescription  for  the  creation  of  holes  during  optimisation,  even  if  these 
holes  would  reduce  the  objective  function.  For  this  reason,  we  believe  that  the  invariant  measure  of  the  method 
described  here  is  of  the  form  (14),  with  p°(f2)  =  0  if  is  not  simply  connected.  Generalisation  of  the  method  to 
include  shapes  with  holes  is  an  important  direction  for  future  work. 


V.  CONCLUSIONS 


We  have  introduced  a  stochastic  level-set  shape  optimisation  method,  which  is  based  on  the  deterministic  (steepest- 
descent)  method  of  [7j  ;I5].  The  stochastic  element  of  the  algorithm  acts  on  the  boundary  of  the  shape  f2,  and  the 
method  converges  to  a  steady  state  in  which  it  explores  a  range  of  shapes,  according  to  a  probability  distribution 
(14).  The  method  is  novel  -  we  are  not  able  to  prove  rigorously  that  it  converges  to  (14)  but  we  have  verified  that 
this  convergence  does  hold  in  two  simple  problems:  shape  matching  (Fig.  [4])  and  perimeter  minimisation  at  fixed  area 
(Fig#  A  deeper  mathematical  analysis  of  the  method  would  be  an  interesting  direction  for  future  work. 

The  motivation  for  introducing  the  stochastic  method  is  to  enable  optimisation  in  non-convex  problems,  in  which 
deterministic  methods  converge  to  local  minima  but  are  not  able  to  search  the  whole  parameter  space  in  order  to  find 
the  global  optimum.  To  demonstrate  the  idea,  we  have  shown  results  in  Fig.  [7]  for  a  simple  non-convex  problem,  in 
which  the  stochastic  algorithm  outperforms  the  deterministic  one  and  converges  to  the  local  optimum.  For  complex 
problems  with  multiple  optima,  an  important  feature  of  the  method  is  that  the  underlying  Boltzmann  distribution  in 
( 14 )  allows  it  to  be  combined  with  methods  such  as  parallel  tempering,  which  are  known  to  be  effective  in  highly  non- 
convex  problems  9,  16,  21.  Finally,  we  considered  a  model  engineering  problem  (Fig.  [8]),  for  which  we  demonstrated 
that  optimisation  can  be  performed  even  for  problems  with  complicated  objective  functions,  although  numerical  errors 
are  signficant  in  that  case. 

Compared  with  other  stochastic  optimisation  methods  that  do  not  use  gradient  (or  sensitivity)  information  (l7j  1221 
173],  the  scheme  presented  here  has  two  strengths.  First,  as  the  noise  strength  is  reduced,  it  recovers  to  the  standard 
deterministic  method,  so  it  is  guaranteed  to  peform  no  worse  than  that  method  (this  property  is  not  guaranteed  in 
many  stochastic  algorithms).  Second,  the  fact  that  the  invariant  measure  (or  at  least  its  T-dependence)  is  known  to 
be  ( 14 1  allows  integration  with  parallel  tempering  and  other  methods  from  statistical  physics [5j  ilBl  1211.  which  have 
proven  extremely  valuable  in  that  context. 
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Appendix  A:  Sensitivity  for  the  curve  perimeter,  local  curvature  and  level-set  reinitialisation 


1.  Local  Curvature 


Several  applications  within  this  work  make  use  of  the  sensitivity  function  associated  with  the  perimeter  of  a  curve. 
The  setting  is  illustrated  in  Fig.  [9j  Recall  the  sensitivity  is  defined  by  ([3|.  For  a  curve  that  is  represented  by  a  set  of 
discrete  boundary  points,  can  be  estimated  as  shown  in  Fig.  [9}  one  displaces  a  boundary  point  a  by  a  distance  S  in 
the  direction  normal  to  the  curve,  and  calculates  the  change  in  F,  which  in  this  case  is  the  perimeter.  It  is  convenient 
also  to  displace  the  boundary  point  by  —  <5,  which  yields  a  central  difference  estimate  of  the  sensitivity:  for  a  boundary 
point  with  index  B ,  this  estimate  is 

sfb  =  ^[F(X1,...,XB+na6,...,Xn)-F(X1,...,XB  -  na6, . . .  Xn)]  (Al) 


where  £B  is  the  length  of  the  boundary  segment  associated  with  point  B1  and  F{X\, . . . ,  Xn)  is  the  value  of  the 
objective  function  given  the  boundary  point  positions  (Xi, . . . ,  Xn). 

When  the  function  F  is  the  perimeter  of  the  curve,  a  short  calculation  shows  that  sB  defined  in  this  way  is  equal 
to  kb  =  ±1  /R,  where  R  is  the  radius  of  the  circle  shown  in  Fig.  [9]^a) .  This  radius  can  be  evaluated  as  shown  in  that 
Figure.  The  sign  of  kb  depends  on  the  direction  of  the  outward  normal  n  at  the  point  B :  for  a  (locally)  convex  shape 
then  nB  >  0  while  for  concave  shapes  kb  <  0. 

It  follows  that  for  deterministic  minimisation  of  shape  perimeter  in  the  absence  of  any  constraint,  one  should  recover 
curvature-driven  flow  of  the  boundary  T  (effects  of  stochastic  noise  on  this  process  have  been  considered  by  Souganidis 
and  Yip  EDI)-  In  our  scheme,  the  computational  implementation  of  this  process  is  affected  by  the  discretisation  of 
the  level-set  field,  and  use  of  discrete  boundary  points.  In  particular,  this  discretisation  affects  our  estimation  of 
sensitivities  [|j.  In  Fig.  |9j(c),  we  initialise  the  level-set  field  as  the  signed  distance  from  a  circle  of  radius  R.  We  then 
infer  the  boundary  points  (as  described  in  Sec.  II B)  and  we  measure  the  curvature  na  at  each  such  point,  using  the 
two  methods  shown  in  Fig.  [8])a,b) .  The  plot  shows  the  average  of  the  Ka  over  the  boundary  points  a.  This  result  is 
compared  with  the  inverse  radius  of  the  circle  1/R.  The  agreement  is  good,  showing  that  the  method  is  valid, 


2.  Level-set  reinitialisation  and  its  effects  on  curvature  measurements 


However,  in  Fig.  [lO]  we  investigate  two  other  effects,  which  lead  to  uncertainties  in  the  curvature.  First,  note 
from  Sec.  II  D|  that  the  level-set  variables  are  periodically  reinitialised,  to  maintain  the  property  that  <j>i  is  the  signed 
distance  of  node  i  from  the  boundary  F.  The  reinitialisation  is  based  on  the  set  of  boundary  point  positions  {Xa}. 
However,  if  one  startes  with  a  set  of  boundary  points,  reinitialises  the  level  set,  and  then  recalculates  a  new  set  of 
boundary  point  positions  {X^},  then  one  does  not  have  in  general  Xa  =  X'a:  that  is,  the  boundary  points  are 
changed  by  the  reinitialisation. 

This  motivates  the  following  numerical  experiment,  whose  results  are  shown  in  Fig.  |10[a).  Start  with  a  level  set 
that  encodes  perfectly  the  signed  distance  from  a  circle  of  radius  R.  Calculate  the  boundary  points  Xa  and  their 
associated  local  curvatures  na.  Then  reinitialise  the  level  set  based  on  the  boundary  points  Xa.  Recalculate  the  new 
set  of  boundary  points  X'a  and  their  associated  curvatures  k'„.  In  the  absence  of  discretisation  errors,  one  would  have 
Ka  =  K'a  =  1  / R  for  all  a.  Fig.  [l0|a)  shows  the  average  of  the  Ka  and  the  Kra,  evaluated  by  summing  over  boundary 
points  and  dividing  by  the  number  of  boundary  points.  One  sees  that  the  average  na  is  close  to  1/R  (as  expected) 
but  the  averaged  K'a  already  shows  significant  numerical  errors  (deviations  from  1  /R). 

Moreover,  Fig.  ©b)  shows  that  the  local  curvatures  evaluated  for  individual  boundary  points  show  significant 
(and  systematic)  variation  according  to  their  positions  on  the  perimeter  of  the  circle.  This  effect  can  be  seen  even 
when  the  level-set  function  is  exactly  equal  to  the  signed  distance  from  a  circle:  for  different  boundary  points,  the 
estimated  sensitivities  are  different  from  1  /R,  indicating  that  inferring  the  discrete  boundary  points  from  the  level  set 
variables  can  affect  the  local  curvature.  This  effect  is  strongest  when  the  boundary  intersects  the  grid  axes  at  angles 
of  15  —  20  degrees.  Moreover,  after  reinitialisation  of  the  level  set  field,  a  similar  effect  is  observed,  but  the  magnitude 
is  strongly  enhaned.  In  this  case  the  effect  is  most  severe  where  the  boundary  T  is  intersecting  the  grid  axes  at  angles 
of  approximately  45  degrees.  The  range  of  k  is  large  enough  in  this  case  that  some  values  are  even  negative. 


3.  Interpretation 

There  are  two  conclusions  from  this  analysis.  First,  even  if  the  level-set  function  is  an  accurate  description  of  a 
shape  (as  with  the  analytic  signed  distance  function  considered  here),  one  expects  uncertainties  in  sensitivities,  due 
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to  discretisation  errors.  This  can  effect  the  convergence  of  the  deterministic  method  to  a  minimum  of  F,  and  the 
extent  to  which  the  stochastic  method  samples  ( 14 ) .  Accurate  estimation  of  sensitivities  is  therefore  an  important 
part  of  any  future  application  of  this  method.  Second,  reinitialisation  of  the  level  set  can  lead  to  small  but  significant 
movements  of  boundary  points,  which  are  large  enough  that  local  sensitivities  change  considerably.  In  determinstic 
optimisation,  the  frequency  of  reinitialisation  reduces  as  the  system  converges  to  the  optimum,  and  this  effect  is  not 
too  pronounced.  On  the  other  hand,  in  stochastic  optimisation,  reinitialisation  is  more  frequent,  and  can  affect  the 
shapes  generated  by  the  method.  A  more  detailed  analysis  of  these  effects  is  a  direction  for  future  work. 
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